What is Homography

单应(Homography)是射影几何中的概念，又称为射影变换。它把一个射影平面上的点(三维齐次矢量)映射到另一个射影平面上，并且把直线映射为直线，具有保线性质。总的来说，单应是关于三维齐次矢量的一种线性变换，可以用一个3×3的非奇异矩阵H表示

Consider two images of a plane (top of the book) shown in Figure 1. The red dot represents the same physical point in the two images. In computer vision jargon we call these corresponding points. Figure 1. shows four corresponding points in four different colors — red, green, yellow and orange. A Homography is a transformation ( a 3×3 matrix ) that maps the points in one image to the corresponding points in the other image.

$\begin{bmatrix} x_2 \\ y_2 \\ 1 \end{bmatrix} = H \begin{bmatrix} x_1 \\ y_1 \\ 1 \end{bmatrix} = \begin{bmatrix} h_{00} & h_{01} & h_{02} \\ h_{10} & h_{11} & h_{12} \\ h_{20} & h_{21} & h_{22} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ y_1 \\ 1 \end{bmatrix}$

这是一个齐次坐标的等式，H乘以一个非零的比例因子上述等式仍然成立，即H是一个3×3齐次矩阵，具有8个未知量.

单应矩阵求解

$\left(\begin{array}{c}x_2\\y_2\\1\end{array}\right)=\left(\begin{array}{ccc}H_{11}&H_{12}&H_{13}\\H_{21}&H_{22}&H_{23}\\H_{31}&H_{32}&H_{33}\end{array}\right)\left(\begin{array}{c}x_1\\y_1\\1\end{array}\right)\Leftrightarrow p_2= Hp_1$

将矩阵的乘法展开，即可得到

$\left\{ \begin{array}{c} x_2 = H_{11}x_1 + H_{12}y_1 + H_{13} \\ y_2 = H_{21}x_1 + H_{22}y_1 + H_{23} \\ 1 = H_{31}x_1 + H_{32}y_1 + H_{33} \end{array} \right.$

方便求解，可以将上面等式变换为Ax=b的形式，做如下变换,第一和第二个式子的左右两边同时乘以第三个式子的左右两边得到

$x_2(H_{31}x_1 + H_{32}y_1 + H_{33})=H_{11}x_1 + H_{12}y_1 + H_{13} \\ y_2(H_{31}x_1 + H_{32}y_1 + H_{33})=H_{21}x_1 + H_{22}y_1 + H_{23}$

将式子的右边变为0

$x_2(H_{31}x_1 + H_{32}y_1 + H_{33})-H_{11}x_1 + H_{12}y_1 + H_{13} = 0 \\ y_2(H_{31}x_1 + H_{32}y_1 + H_{33})-H_{21}x_1 + H_{22}y_1 + H_{23} = 0$

将上面的等式改写为向量积的形式，令

$h=(H_{11}, H_{12},H_{13},H_{21},H_{22},H_{23},H_{31},H_{32},1)^T$

单应矩阵H是一个齐次矩阵，可以将其最后一个元素归一化为1。

$H_{11} x_1 + h_{12} y_1 + H_{13} - H_{31} x_1 x_2 - H_{32}y_1 x_2 = x_2$ $H_{21} x_1 + h_{22} y_1 + H_{23} - H_{31} x_1 y_2 - H_{32}y_1 y_2 = y_2$

一对匹配的点对，可以得到上述等式，H有8个未知量，也就说最少4对匹配的点对(任意3点不共线)，就可以求出两幅图像的单应矩阵H。但是通常来说，图像的匹配点对要超过4对，设得到了n对匹配的点对，可以得到如下的等式

对于大于4对点的情况可以使用线性最小二乘求解。

一个简单的例子

/**
******************************************************************************
* @file
* @author  maky <chengwei920412@outlook.com>
* @version
* @date    2018-01-04 16:22:58
* @brief
******************************************************************************
* @attention
*
*
******************************************************************************
*/
#include <iostream>
#include <Eigen/Eigen>
#include <Eigen/geometry>
#include <opencv2/opencv.hpp>
#include <opencv2/aruco.hpp>
#include <opencv2/core/eigen.hpp>

cv::Mat findHomography(const std::vector<cv::Point2f> &xs1, const std::vector<cv::Point2f> &xs2)
{
    if (xs1.size() != xs2.size()) {
        return cv::Mat();
    }
    int size = xs1.size();
    Eigen::Matrix<double, Eigen::Dynamic, 9> A = Eigen::MatrixXd::Zero(2 * size, 9);
    for (auto pos = 0; pos < size; pos++) {
        auto &x1 = xs1.at(pos);
        auto &x2 = xs2.at(pos);
        auto rowx = pos * 2;
        auto rowy = rowx + 1;
        A(rowx, 0) = x1.x;
        A(rowx, 1) = x1.y;
        A(rowx, 2) = 1;
        A(rowx, 6) = -1 * x1.x * x2.x;
        A(rowx, 7) = -1 * x1.y * x2.x;

        A(rowy, 3) = x1.x;
        A(rowy, 4) = x1.y;
        A(rowy, 5) = 1;
        A(rowy, 6) = -1 * x1.x * x2.y;
        A(rowy, 7) = -1 * x1.y * x2.y;

        A(rowx, 8) = -1 * x2.x;
        A(rowy, 8) = -1 * x2.y;
    }
    Eigen::JacobiSVD<Eigen::Matrix<double, Eigen::Dynamic, 9>> solver(A, Eigen::ComputeFullV);
    const Eigen::VectorXd x = solver.matrixV().rightCols(1);
    Eigen::Matrix3d h;
    h << x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8];
    h.normalize();
    cv::Mat homography;
    cv::eigen2cv(h, homography);
    return homography;
}

int main(int argc, const char **argv)
{
    cv::Mat source = cv::imread("1546583988066087246-1.png");
    cv::Mat patch = cv::imread("patch.png");
    auto dictionary = cv::aruco::getPredefinedDictionary(cv::aruco::DICT_6X6_250);
    cv::Ptr<cv::aruco::DetectorParameters> parameters = cv::aruco::DetectorParameters::create();
    // detect
    std::vector<int> ids;
    std::vector<std::vector<cv::Point2f>> corners, rejected_candidates;
    cv::aruco::detectMarkers(source, dictionary, corners, ids, parameters, rejected_candidates);

    std::vector<cv::Point2f> patch_points;
    patch_points.push_back(cv::Point2f(0., 0.));
    patch_points.push_back(cv::Point2f(patch.cols - 1, 0.));
    patch_points.push_back(cv::Point2f(patch.cols - 1, patch.rows - 1));
    patch_points.push_back(cv::Point2f(0., patch.rows - 1));

    //cv::Mat homography = cv::findHomography(patch_points, corners[0]);
    cv::Mat homography = findHomography(patch_points, corners[0]);
    cv::Mat output;
    cv::warpPerspective(patch, output, homography, source.size());

    for (auto u = 0; u < source.cols; u++) {
        for (auto v = 0; v < source.rows; v++) {
            auto pixel = output.at<cv::Vec3b>(v, u);
            if (pixel[0] != 0 || pixel[1] != 0 || pixel[2] != 0) {
                source.at<cv::Vec3b>(v, u) = pixel;
            }
        }
    }
    cv::imwrite("patched.png", source);
    cv::imshow("source", source);
    cv::imshow("patch", patch);
    cv::imshow("output", output);

    while (true) {
        auto key = cv::waitKey(1);
        if (' ' == key) {
            break;
        }
    }
    return 0;
}